% Author:   Vasyl Mykhalchuk
% Date:     22/04/2012

% v_initial, v_deformed - row vectors, i.e. dimensions [1, 3]
%{
v_initial - initial coordinate of the vertex
v_deformed - coordinate of the vertex after the deformation
neighbor_list - vertex adjacency matrix (CCW)

%}
function [ strains ] = strain_vert( v_initail, v_deformed, neighbor_list, frame_initial, frame_deformed )
disp = v_deformed - v_initail;          % Point displacement
strains = zeros(1, 2);                  % Strain row vector: [ average_principal_strain_max, average_principal_strain_min]

neighb_num = size(neighbor_list, 2);
nonzero_counter = 0;
for i = 1 : neighb_num
    if (neighbor_list(i) ~= 0)
        nonzero_counter = nonzero_counter + 1;    
    else
        i = neighb_num + 1;
        break;
    end % IF
end % FOR

neighb_num = nonzero_counter;

ae = zeros(neighb_num, 3);              % Adjacent edge vectors initial
ae_p = zeros(neighb_num, 3);            % Adjacent edge vectors deformed

% Compute adjacent edge vectors
for i = 1 : neighb_num
    adjv_id = neighbor_list(i);
    
    ae(i, :) = frame_initial(adjv_id, :) - v_initail; 
    ae_p(i, :) = frame_deformed(adjv_id, :) - v_deformed;
    
end % FOR

% Strains measured from edge length differencee
e_p = zeros(neighb_num, 1);
e_q = zeros(neighb_num, 1);
e_r = zeros(neighb_num, 1);

% Angles between the 1st and the 2nd(q), the 1st and the 3rd(r) elements in
% the strain gage
ang_q = zeros(neighb_num, 1);
ang_r = zeros(neighb_num, 1);

% Normal and Shear strains computed for each pair of adjacent triangles
e_x = zeros(neighb_num, 1);
e_y = zeros(neighb_num, 1);
e_xy = zeros(neighb_num, 1);
phi = zeros(neighb_num, 1); 

% Normal vector for each adjacent edge vector
n = zeros(neighb_num, 3);

e1 = zeros(neighb_num, 3);
e2 = zeros(neighb_num, 3);
e1_m = zeros(1, 3);
e2_m = zeros(1, 3);

%
for j = 1 : neighb_num
    % 1. Compute eps_x, eps_y, eps_xy for each adjacent triangle
    % A, B, C - indices of the 3 elements of the strain gage
    A = j;
    B = mod( j + 1, neighb_num + 1);
    C = mod( j + 2, neighb_num + 1);
    if (B == 0)
        B = 1;
        C = 2;
    elseif (C == 0)
        C = 1;
    end % IF
    
    e_p(j) = norm( ae_p(A, :) - ae(A, :) ) / norm( ae(A, :) );
    e_q(j) = norm( ae_p(B, :) - ae(B, :) ) / norm( ae(B, :) );
    e_r(j) = norm( ae_p(C, :) - ae(C, :) ) / norm( ae(C, :) );
    
    % Normalized 3 edge vectors
    ae_a = ae_p(A, :) / norm( ae_p(A, :) );
    ae_b = ae_p(B, :) / norm( ae_p(B, :) );
    ae_c = ae_p(C, :) / norm( ae_p(C, :) );
    
    cos_q = dot(ae_a, ae_b);
    cos_r = dot(ae_a, ae_c);
    % Angles between the edges
    ang_q(j) = acos(cos_q);
    ang_r(j) = acos(cos_r);
    
    % Compute principal strains and the orientation
    left = zeros(2, 1); right = zeros(2, 2); unknowns = zeros(2, 1);
    left(1, 1) = e_q(j) - 0.5 * (1 + cos(2 * ang_q(j)) ) * e_p(j);
    left(2, 1) = e_r(j) - 0.5 * (1 + cos(2 * ang_r(j)) ) * e_p(j);
    right(1, 1) = 0.5 * (1 - cos(2 * ang_q(j)) );
    right(1, 2) = 0.5 * sin(2 * ang_q(j));
    right(2, 1) = 0.5 * (1 - cos(2 * ang_r(j)) );
    right(2, 2) = 0.5 * sin(2 * ang_r(j));
    unknowns = right \ left;
    e_y(j) = unknowns(1, 1);
    e_xy(j) = unknowns(2, 1);
    
    % Now the normal strain and shear strain are computed for the
    % rectangular coordinate system
    % Let's now compute the principal strains e1 and e2 from e_x, e_y, e_xy
    % First, compute e_1, e_2, and phi (all scalars) using Mohr's circle
    R = sqrt( e_xy(j) ^ 2 + (e_x(j) - e_y(j)) ^ 2 ) / 2;
    e_1(j) = ( e_x(j) + e_y(j) ) / 2 + R;
    e_2(j) = ( e_x(j) + e_y(j) ) / 2 - R;
    phi(j) = atan( e_xy(j) / (e_x(j) - e_y(j)) ) / 2;
    
    % Having computed phi, we can compute e1 by rotating e_x by phi
    n(j, :) = cross(ae_a, ae_b);
    % Rotate ae_p vector around n by the angle phi
    % r' = (costh) r + (1 - costh) n (n.r) + (sinth) n x r
    cphi = cos(phi(j));
    sphi = sin(phi(j));
    n_dot_r = dot(ae_p(j, :), n(j, :));
    n_crs_r = cross(ae_p(j, :), n(j, :));
    e1(j, :) = cphi * ae_p(j, :) + (1 - cphi) * n_dot_r * n(j, :) + sphi * n_crs_r;
    % e2 is computed by rotating e1
    e2(j, :) = cross(n(j, :), e1(j, :));
end % FOR

% Average principal strains and orientation
eps1 = mean(e_1);
eps2 = mean(e_2);

strains = [eps1, eps2];

end % FUNCTION strain_vert

